{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GPyTorch Regression Tutorial\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this notebook, we demonstrate many of the design features of GPyTorch using the simplest example, training an RBF kernel Gaussian process on a simple function. We'll be modeling the function\n",
    "\n",
    "\\begin{align}\n",
    "y &= \\sin(2\\pi x) + \\epsilon \\\\\n",
    "  \\epsilon &\\sim \\mathcal{N}(0, 0.04) \n",
    "\\end{align}\n",
    "\n",
    "with 100 training examples, and testing on 51 test examples.\n",
    "\n",
    "**Note:** this notebook is not necessarily intended to teach the mathematical background of Gaussian processes, but rather how to train a simple one and make predictions in GPyTorch. For a mathematical treatment, Chapter 2 of Gaussian Processes for Machine Learning provides a very thorough introduction to GP regression (this entire text is highly recommended): http://www.gaussianprocess.org/gpml/chapters/RW2.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up training data\n",
    "\n",
    "In the next cell, we set up the training data for this example. We'll be using 100 regularly spaced points on [0,1] which we evaluate the function on and add Gaussian noise to get the training labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training data is 100 points in [0,1] inclusive regularly spaced\n",
    "train_x = torch.linspace(0, 1, 100)\n",
    "# True function is sin(2*pi*x) with Gaussian noise\n",
    "train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * math.sqrt(0.04)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the model\n",
    "\n",
    "The next cell demonstrates the most critical features of a user-defined Gaussian process model in GPyTorch. Building a GP model in GPyTorch is different in a number of ways.\n",
    "\n",
    "First in contrast to many existing GP packages, we do not provide full GP models for the user. Rather, we provide *the tools necessary to quickly construct one*. This is because we believe, analogous to building a neural network in standard PyTorch, it is important to have the flexibility to include whatever components are necessary. As can be seen in more complicated examples, this allows the user great flexibility in designing custom models.\n",
    "\n",
    "For most GP regression models, you will need to construct the following GPyTorch objects:\n",
    "\n",
    "1. A **GP Model** (`gpytorch.models.ExactGP`) -  This handles most of the inference.\n",
    "1. A **Likelihood** (`gpytorch.likelihoods.GaussianLikelihood`) - This is the most common likelihood used for GP regression.\n",
    "1. A **Mean** - This defines the prior mean of the GP.(If you don't know which mean to use, a `gpytorch.means.ConstantMean()` is a good place to start.)\n",
    "1. A **Kernel** - This defines the prior covariance of the GP.(If you don't know which kernel to use, a `gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())` is a good place to start).\n",
    "1. A **MultivariateNormal** Distribution (`gpytorch.distributions.MultivariateNormal`) - This is the object used to represent multivariate normal distributions.\n",
    "  \n",
    "  \n",
    "### The GP Model\n",
    "  \n",
    "The components of a user built (Exact, i.e. non-variational) GP model in GPyTorch are, broadly speaking:\n",
    "\n",
    "1. An `__init__` method that takes the training data and a likelihood, and constructs whatever objects are necessary for the model's `forward` method. This will most commonly include things like a mean module and a kernel module.\n",
    "\n",
    "2. A `forward` method that takes in some $n \\times d$ data `x` and returns a `MultivariateNormal` with the *prior* mean and covariance evaluated at `x`. In other words, we return the vector $\\mu(x)$ and the $n \\times n$ matrix $K_{xx}$ representing the prior mean and covariance matrix of the GP. \n",
    "\n",
    "This specification leaves a large amount of flexibility when defining a model. For example, to compose two kernels via addition, you can either add the kernel modules directly:\n",
    "\n",
    "```python\n",
    "self.covar_module = ScaleKernel(RBFKernel() + WhiteNoiseKernel())\n",
    "```\n",
    "\n",
    "Or you can add the outputs of the kernel in the forward method:\n",
    "\n",
    "```python\n",
    "covar_x = self.rbf_kernel_module(x) + self.white_noise_module(x)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will use the simplest form of GP model, exact inference\n",
    "class ExactGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
    "    \n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "# initialize likelihood and model\n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = ExactGPModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model modes\n",
    "\n",
    "Like most PyTorch modules, the `ExactGP` has a `.train()` and `.eval()` mode.\n",
    "- `.train()` mode is for optimizing model hyperameters.\n",
    "- `.eval()` mode is for computing predictions through the model posterior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training the model\n",
    "\n",
    "In the next cell, we handle using Type-II MLE to train the hyperparameters of the Gaussian process.\n",
    "\n",
    "The most obvious difference here compared to many other GP implementations is that, as in standard PyTorch, the core training loop is written by the user. In GPyTorch, we make use of the standard PyTorch optimizers as from `torch.optim`, and all trainable parameters of the model should be of type `torch.nn.Parameter`. Because GP models directly extend `torch.nn.Module`, calls to methods like `model.parameters()` or `model.named_parameters()` function as you might expect coming from PyTorch.\n",
    "\n",
    "In most cases, the boilerplate code below will work well. It has the same basic components as the standard PyTorch training loop:\n",
    "\n",
    "1. Zero all parameter gradients\n",
    "2. Call the model and compute the loss\n",
    "3. Call backward on the loss to fill in gradients\n",
    "4. Take a step on the optimizer\n",
    "\n",
    "However, defining custom training loops allows for greater flexibility. For example, it is easy to save the parameters at each step of training, or use different learning rates for different parameters (which may be useful in deep kernel learning for example)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/50 - Loss: 0.939   lengthscale: 0.693   noise: 0.693\n",
      "Iter 2/50 - Loss: 0.908   lengthscale: 0.644   noise: 0.644\n",
      "Iter 3/50 - Loss: 0.874   lengthscale: 0.598   noise: 0.598\n",
      "Iter 4/50 - Loss: 0.837   lengthscale: 0.555   noise: 0.554\n",
      "Iter 5/50 - Loss: 0.795   lengthscale: 0.514   noise: 0.513\n",
      "Iter 6/50 - Loss: 0.749   lengthscale: 0.476   noise: 0.474\n",
      "Iter 7/50 - Loss: 0.699   lengthscale: 0.440   noise: 0.437\n",
      "Iter 8/50 - Loss: 0.649   lengthscale: 0.405   noise: 0.402\n",
      "Iter 9/50 - Loss: 0.600   lengthscale: 0.372   noise: 0.369\n",
      "Iter 10/50 - Loss: 0.556   lengthscale: 0.342   noise: 0.339\n",
      "Iter 11/50 - Loss: 0.516   lengthscale: 0.315   noise: 0.310\n",
      "Iter 12/50 - Loss: 0.480   lengthscale: 0.291   noise: 0.284\n",
      "Iter 13/50 - Loss: 0.448   lengthscale: 0.270   noise: 0.259\n",
      "Iter 14/50 - Loss: 0.413   lengthscale: 0.254   noise: 0.237\n",
      "Iter 15/50 - Loss: 0.380   lengthscale: 0.241   noise: 0.216\n",
      "Iter 16/50 - Loss: 0.355   lengthscale: 0.231   noise: 0.197\n",
      "Iter 17/50 - Loss: 0.314   lengthscale: 0.223   noise: 0.179\n",
      "Iter 18/50 - Loss: 0.292   lengthscale: 0.218   noise: 0.163\n",
      "Iter 19/50 - Loss: 0.262   lengthscale: 0.214   noise: 0.148\n",
      "Iter 20/50 - Loss: 0.236   lengthscale: 0.214   noise: 0.135\n",
      "Iter 21/50 - Loss: 0.201   lengthscale: 0.216   noise: 0.122\n",
      "Iter 22/50 - Loss: 0.176   lengthscale: 0.220   noise: 0.111\n",
      "Iter 23/50 - Loss: 0.158   lengthscale: 0.224   noise: 0.102\n",
      "Iter 24/50 - Loss: 0.125   lengthscale: 0.231   noise: 0.093\n",
      "Iter 25/50 - Loss: 0.101   lengthscale: 0.239   noise: 0.085\n",
      "Iter 26/50 - Loss: 0.078   lengthscale: 0.247   noise: 0.077\n",
      "Iter 27/50 - Loss: 0.066   lengthscale: 0.256   noise: 0.071\n",
      "Iter 28/50 - Loss: 0.052   lengthscale: 0.265   noise: 0.065\n",
      "Iter 29/50 - Loss: 0.036   lengthscale: 0.276   noise: 0.060\n",
      "Iter 30/50 - Loss: 0.036   lengthscale: 0.286   noise: 0.056\n",
      "Iter 31/50 - Loss: 0.031   lengthscale: 0.297   noise: 0.052\n",
      "Iter 32/50 - Loss: 0.028   lengthscale: 0.306   noise: 0.048\n",
      "Iter 33/50 - Loss: 0.030   lengthscale: 0.315   noise: 0.045\n",
      "Iter 34/50 - Loss: 0.035   lengthscale: 0.322   noise: 0.043\n",
      "Iter 35/50 - Loss: 0.039   lengthscale: 0.326   noise: 0.041\n",
      "Iter 36/50 - Loss: 0.043   lengthscale: 0.329   noise: 0.039\n",
      "Iter 37/50 - Loss: 0.047   lengthscale: 0.327   noise: 0.038\n",
      "Iter 38/50 - Loss: 0.052   lengthscale: 0.323   noise: 0.037\n",
      "Iter 39/50 - Loss: 0.048   lengthscale: 0.317   noise: 0.036\n",
      "Iter 40/50 - Loss: 0.051   lengthscale: 0.309   noise: 0.036\n",
      "Iter 41/50 - Loss: 0.051   lengthscale: 0.302   noise: 0.036\n",
      "Iter 42/50 - Loss: 0.047   lengthscale: 0.295   noise: 0.036\n",
      "Iter 43/50 - Loss: 0.048   lengthscale: 0.288   noise: 0.036\n",
      "Iter 44/50 - Loss: 0.047   lengthscale: 0.281   noise: 0.037\n",
      "Iter 45/50 - Loss: 0.047   lengthscale: 0.276   noise: 0.037\n",
      "Iter 46/50 - Loss: 0.040   lengthscale: 0.273   noise: 0.038\n",
      "Iter 47/50 - Loss: 0.037   lengthscale: 0.271   noise: 0.039\n",
      "Iter 48/50 - Loss: 0.040   lengthscale: 0.270   noise: 0.040\n",
      "Iter 49/50 - Loss: 0.033   lengthscale: 0.269   noise: 0.042\n",
      "Iter 50/50 - Loss: 0.032   lengthscale: 0.269   noise: 0.043\n"
     ]
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "training_iter = 2 if smoke_test else 50\n",
    "\n",
    "\n",
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
    "], lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "for i in range(training_iter):\n",
    "    # Zero gradients from previous iteration\n",
    "    optimizer.zero_grad()\n",
    "    # Output from model\n",
    "    output = model(train_x)\n",
    "    # Calc loss and backprop gradients\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/%d - Loss: %.3f   lengthscale: %.3f   noise: %.3f' % (\n",
    "        i + 1, training_iter, loss.item(),\n",
    "        model.covar_module.base_kernel.lengthscale.item(),\n",
    "        model.likelihood.noise.item()\n",
    "    ))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make predictions with the model\n",
    "\n",
    "In the next cell, we make predictions with the model. To do this, we simply put the model and likelihood in eval mode, and call both modules on the test data.\n",
    "\n",
    "Just as a user defined GP model returns a `MultivariateNormal` containing the prior mean and covariance from forward, a trained GP model in eval mode returns a `MultivariateNormal` containing the posterior mean and covariance. Thus, getting the predictive mean and variance, and then sampling functions from the GP at the given test points could be accomplished with calls like:\n",
    "\n",
    "```python\n",
    "f_preds = model(test_x)\n",
    "y_preds = likelihood(model(test_x))\n",
    "\n",
    "f_mean = f_preds.mean\n",
    "f_var = f_preds.variance\n",
    "f_covar = f_preds.covariance_matrix\n",
    "f_samples = f_preds.sample(sample_shape=torch.Size(1000,))\n",
    "```\n",
    "\n",
    "The `gpytorch.settings.fast_pred_var` context is not needed, but here we are giving a preview of using one of our cool features, getting faster predictive distributions using [LOVE](https://arxiv.org/abs/1803.06058)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get into evaluation (predictive posterior) mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Test points are regularly spaced along [0,1]\n",
    "# Make predictions by feeding model through likelihood\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    test_x = torch.linspace(0, 1, 51)\n",
    "    observed_pred = likelihood(model(test_x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the model fit\n",
    "\n",
    "In the next cell, we plot the mean and confidence region of the Gaussian process model. The `confidence_region` method is a helper method that returns 2 standard deviations above and below the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAADDCAYAAABtec/IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXl4U2X2+D/ZmnRPW2gLZWtBAUFtaVDBDWhRAdcKwvh1dGYQdAac34wOKqOizig4IoM7DoI6DuJWRR3BhUVFBISUFsGhsrSUrXub7k2z3N8fadI0TdrQpm1K3s/z8JDe9eTm3nPPe96zyCRJQiAQBCby3hZAIBD0HkIBCAQBjFAAAkEAIxSAQBDACAUgEAQwyq4eQKfTpTd/nKrX6x/q6vEEAkHP0SULoPnhn6XX67cA43Q63TjfiCUQCHoCma/iAHQ63TG9Xj/cJwcTCAQ9gk98ADqd7kHgHl8cSyAQ9By+tAA+BObp9XqDu/UPP/ywCDkUCHqJZ555RuZueZecgPYxv16v3wfkAfOBZz1t/+STT3Z4zJKSEmJjY7siVrfj7zL6u3zg/zL6u3zgvYyPP/64x3VdHQKkA9HNn7XYlIBAIOgjdFUBrAaSdDrdfAC9Xp/ZdZEEAkFP0aUhQPN4f7WPZBEECGazmZqaGmpqavDXbFSr1Up1dXVvi9EurjLKZDLUajXx8fEold492l0OBBIIzpaioiIiIyOJiYlBJnPrm+p1TCYTKpWqt8VoF1cZJUnCYDBQVFTEoEGDvDqGCAUW9DhGo5GIiIhef/izs7PJzs7u9vMYDAY+/vjjbj+PTCZDq9ViNBq93kcoAEGPI0mS1w9/YWEh6enpFBUVdfp82dnZrFmzhq1bt7JmzRry8my+6sjISDIzu99tpdVq3Z4nOzub0aNH8/HHH/Pxxx+zYsUKh2zuaG+dHZlMdlbDKjEEEPg1y5YtY+fOnSxdupQXX3zxrPc3GAwsX76c9evXO5bdfvvtrF+/nujo6Hb29C1RUVFtlqWkpJCYmEhGRoZj2fTp09m0aVObbfPy8li7di1PP/20T+USCkDgl2i1WhobGx1/r169mtWrV6PRaDAY3MaauSUzM5MpU6a0WhYVFcXWrVtJTU0lOzubrVu3kpOTw9y5c8nKygIgKyuLmTNnsm3bNqKjo0lMTCQ/P5/MzEwSExMZOXIkX375JevXr2fBggU88MADAK22T0xMZO3atSQnJ7Nv3z6vv7f9Tb9t2zYApkyZQk5ODvn5+WRnZxMZGcm2bduwWCxMnTqVpKQkr6+HK2IIIPBLDh06xOzZswkODgYgODiYOXPmkJube9bHqqqq8rguJSWFtLQ0kpOTWbt2LTk5OWzbto3Jkyfz6KOPkpqa6nj4p0yZQlRUFE8//TR33nmn4xgZGRkkJSW12f6RRx7hlltuIS0tjcTExLOSOSkpiejoaKKjo9mwYQNTpkwhMTGRlJSUNuu6glAAAr9kwIABREREYDQa0Wg0DsdhfHz8WR1nypQpjre6nfz8fNLS0lotsw8HbrnlFubOncuKFStoamoiMjKSlJQUhxWh1WpbHXvFihWkpqY6lrluf7YYDAaSkpJYsWIFkZGRJCcnO5aDbShgX3fxxRe3WtcZxBBA4LeUlJQwb9485s6dy9q1azvlCExKSmLRokWsWbOGxMREcnJyePnllx3rDQZDqyGA3WSfPHkyU6dOZe3atY63r90ENxgMaLVaZs6cySOPPOJQCk899VSr7R944AE2bNhAcnKyY9+UlBTHubOzs8nPz3fMEOTn5ztks5+vqqqKvLw8KisrMRgM5OfnO9ZVVFSQl5dHfn5+q+OeDT5LBuqIhx9+WBK5AD2Dv8t39OhRhg4d6tfz7H0xDsDO0aNHGTFihOPvxx9/3GMykBgCCAQBjFAAAkEAIxSAQBDACAUgEAQwQgEIBAGMUAACQQAjFIDgnCY7O5sJEya0yvrLy8trsyxQEYFAgl5Fo1H75DiNje5TYFNSUhyBQK+88gpgCw22h9UGOkIBCM55IiMjPa7Ly8trlcDjmmiTn5/PihUreOCBB9i2bZvPs/F6my4PAXQ63fzmf//whUCCwKKx0eiTfx2RkZHBmjVr2oTjuibwuCbapKWlodVqSUtL61LMvb/ii9ZgW/R6vb04aHpH+wgEvUFaWpojvdYV5wQed4k27nL5zxW6agEkYSsNDraS4J1PTBYIuoHs7GzWrl1LXl6e401vLwWWnZ3tSODZunUrFRUVDktg//795OXlsWnTJvLz8x1JN+ea47CrVYGdKwKPA97vmjgCgW9JSUlxVAOyF+1ISUnh0KFDjm2cx/X24homk4lZs2YBtgpCgNtKPX0dnzgBmzsE7WvuEOSRkpKSDo/VF8ZZ/i6jv8tntVqxWCy9LUa7+Lt84FlGq9Xq1bMGvpsFSNfr9Q91tJG3Kar+nMpqx99l9Gf5qqurUSgUfp9u6+/ygXsZ5XK517+/T2YB9Hr9s82fhRNQIOhD+GIW4B86ne6YTqer9JFMAoGgh+iqE3ALcO7OkQgEHjAYDGzbtq1VSe++iIgEFPQ6L31zrEv73zd5eLvrs7OzycrKclTmPZuIvhUrVpCcnExOTo5jFgFamn0IBSAQ+DHuGoN426bLYDBQUVFBWlqa2yYi50KAkFAAgnOazMzMVmW7wRYPYM8BSExMpKqqisjIyFYx/0888QRZWVnk5+ezdetWHn30UXbt2oXBYGjT7MM1n6CioqJN/oC9s49dluTk5Fb79FZikkgHFgQc9iYed999N2lpaWRmZrqN+bdHDqalpTFu3DgAt80+XPMJ3B1rxYoVzJ07l4yMDKZMmdJmn95CWACCc5qZM2fyhz/8odWyrVu3Ajjq+xsMhi7H/EdGRpKUlORQCu6OZR9G2PsIuO7TGwgFIDin0Wq1rRqDVFVVkZyczFNPPUVmZibR0dFkZGSQn5/fKuY/JyeH6upqRyOQffv2kZ2d7bbZh2tDENdj2fdbvny5463vuo9zx6GeRDQG6QT+LqO/yycag/gG0RhEIBB0CaEABIIARigAQY8jk8noqaFnoCFJEjKZW2vfLcIJ2AGFhYX8+te/Zt26dW1aUxvqTRTXGKluMFFrNFPXZKHOaMZiBbkMZDKQy2RoVAoiNErCNUoig1X0CwsiMti/x5fdiVqtprq6mpiYmLO6WQXtI0kSBoMBtdr7QqtCAXTAsmXL2LlzJ0uXLuWJZc+RV1ZP7olKGg7X02jqfM54mFrJgEgNA7UaEmNCCdcEzk8RHx9PQUEBVVVVfmsJWK1W5HL/NpBdZZTJZKjV6jYvqvYInLvuLNFqtTQ2Njr+Xr16NatXr0apCuLRd3cQGmq7dNUVpaxb9hfuWPwcEdH9vT5+rdHMkZJajpTUsp1y4iPVnB8bxojYMEKCFD7/Pv6EUqkkPDzcr2cq/H0mBXwjo3+ruF7CYpXYuF3PhKk3olJrAFCpNaRMnsHit75ute2W9as4/nMWW9av6vT5JCQKqxr574+H0E28mne+PUBhVWPHOwoEXURYAE40ma3sP1XFgdPV1DXJkFQazE1GlEFqzE1GNCFhRET3p66ujsU3pmA2NTn23bXxfXZtfB+lKohln3WucOSW9avI/zmLNS89R8XCJcRHaBg3REtSvxAxVhZ0C0IB0PLgZ5+swmhuGdfXGiq4bPpsLps+i92bPqSmssyxbvFbX/P568s5uGsrJmMjKrWGsRPTuP7uRWd9/o6USUyomonDoxkWE9K1LyoQuBDQQwCLVWLfCQP/3n2C3fkVrR5+gLsee4GMhY8xMGkUGQsf467HXnCsi4jujzoktI2FAPDqors4fSyXVxfdRXVFaYdyLH7ra1ImzfA43Mg/eYpp117Dm1v3U1ztuQlGYWEh6enpXheEFAgCVgHkl9Wxfs9JfjhWflbe/OqKUt5ccg+nj+Xy0/dfMy7tRu5buZ7Lps+mprLc4RN499kHvfYNeFImdqei/Zhvv/pPPsw6zZZDJTQ0tZXZPmPx/PPPe38hBAGNTxRAc1nwPkFlfROf/VTI5weKMDSYqK4odfum9rR8y/pVnDiUw7vPPkhDbRVBag0Dk0ax9+uPObhzC7s2vo8kSRSfOIYkSeza+D6Lpo1h8Y1t872dz2Efbjgrk8U3prBo2hjHMXdtfJ+/TLuAmy8Zwbo9J8ktqgFsMxYajYbVq1djtVpZt24dGo2m1xJMBH2HLvsAmguD/gtovy5TL2O1SmSdMLC3oBKLtWXu2dmLn7FwicflruP04hO2Mlb28bpcoSQ0MprG+losTtspVEFoQsKYv/T1NjI5n8N5eJGx8DHApiA8+RkaTRY2Hyoht6iGPdkHePqJR/nss89oaGhAo9Fw880388wzz/juAgrOSbqsAPR6/RadTpfnC2G6i+JqI9t+KaWstmX87Mnx5op9uUKpImXSDA7s3IK5qeU4yiA1F16ejlyuIGvrF0A8yNQgmUEGFpOJuqpydm/6wKFgOnL6OccWtDc0ADhZ2UBxtRyLUoPRaESjsf0fERFxVgEhgsDknJ4FMFus/Hi8kuyTbSPOPHnxr771d3yX+Ybbt+7md15tfrtHAONANgFzk47sb4YBCUCc7eD2UzmdctfGYnZtPI5MdoKrbs3l9JHVHD/0GuYmQ5sZBGfLoL2ZCDtNFiu5+ae5JuN2/vrnBbz+2qsUFxf7+GoKzkV6VAH0ZGuw8joT2/OqqGwwu12vUIcgb36rKlVBmJuMKFQatPGD2yyXKzWUnInhSM4MNGHP0FAzBFC0esBtWIBigjQy1MHhNNbXYzKagX7YlEMcknQp32UCLAOeBJkek3EzTU3lLPvNNW0sAwClKohrf3s/E2++k8yVj1B0qoDwqH6tzjzz/qW2fcoU3HHvnxg9uL9fzwb4e/syf5cPfCNjjyqAnmgNZrVK7Dtp4MeCWqxyNaGhrRMjnM1rY211m7draGioY/nI1N/y9bpGftquI2vzQGAsAHKFlYFJDQwd1cDgUQ0c/OEVDu58G4WqEqu5gdS028hYuISPXnqSH7/4ELlSg8UUw5gJv2PsxPmcOhxM1tZKGusTQZoITOTnHyBu6NMoFJmUnHoac9OxVpZBaGgoX725ghOHctj5yb9b+Stc2V1YjzpayaXDovw6gMjfQ239XT7ouozn1BCgqsHE5kMl7YbRduR4s1pg/DVr+eG/0bz1tzDH+jCtmbETqxk7sYa4xFK00S1BOQd/2MuEGVe2MdPbmu87OX/cFPZ89RcWvf4cqiAT+T+HcPCHcA7sjKC4IBz4LfBrZPKPMBmXowkJc2sZtBd1KCGx93glhVWNXDM6llD1OfUzC3xIl0uC6XS6mcDrwDy9Xp/pabvuLgmWW1TDd4fLaLJY3a53dbzZsT9ElSVl/GvxDqzmP1NZEmxbF2QlZVIVuvQqhl1Qj7w5R6euro7Q0NCzlhHg45f/xu5NH3DZ9NtavcXNTTJy9WFseDWf6orLQbKdLCRiPzfeIyd3z2J+3t3WL+EuAanoVAEfv7CEOxY/x4D4AVw7JpYEbXCn5O0u/D3Zxt/lA+9lbK8kmC9mATIBjw9+d2M0Wfj2SBmHi2tbLXfN0vPk9Jv+20Xs/TqST1/rh7HhKgCi45uYeH0F468xEBLuXqGcLR15/pVBEmMn1jB2Yj8qS46x49MYfvxCS331xby3HMKjn8Vk/B3KoBy3swHOfPfh2lZTmJ/kFHLliBguGhTpk+8iOHfo07bhGUMjXx8qoabR1Gad6zy+a7SdyWikrmoKT9+ZgM2DD/A/4Akqij7iy38rufrWziX1uONscgeiYs3cMK+YqbeXsvuLKL75MIaaiguA3SSOOUVo5ApqKg+2PUcHSqak1sik8/qhVARsAKjAhT6pAKxWiT0FlegLDG2n99p5CEaNv4rLps9mRPJcPn0tnsP7bE49dfBpzKZHsZjfRqUOYuzEaZ1K6mmP9sJ9PdUU0IRamTSznEunVfLth/3YviGaI9mDUCj/yeTbyjA1laEKavn+diVjj1VwVTKHCmuorDMxbWwcYcIvIKAP5gJUNZj4KPsMe49Xuq0m4ymxZuHKd6muNGM2LeedZWlUl48lJNzMLQsKSb76L1gt/0YZpOrQvO4s1RWlbnMHoOOaAsGhVqb9poSH1hxFN9WAxSxjy/r+rPxDEsd+anFG2pWMxdTkMXCoqLqRD/Sn200qEgQOfUoB/K+wmvf2nqKo2rOX39Ob9os38zhx6B32fj0ICZh4QwUPrT3KxOsrqasuaxOH72u2rF/VKncgY+Fj5O7d3ibW3zlvwDUfQdvfzOz7z/D75ceJHWyk9LSa1x4axrpnomios/2UtYYKUq/JaPe71DWZ+Tj7DEdKatusEwQWfcIOrG+y8M0vpeSV1Xm1vfP028qFC9i1MQO4rXntTiTrvfz4RS6F+Rdzx+LnWk0Hpt9+L+uW/YXqilKfWAHtDUk68gt4ylPoN7CAkPDbGZi0gDN5t7L/uwEUHIrhV385w12PveCYpbBPbbrDbLXy1c8lVNQ1cWli2863gsDA7y2AY6W2tF1vH36w5fHfsuAxio5fgkJxGNvDXwssRBmUTsrkISRfPc2t2e2LEl/OtJfr78laWfaba9q1DJ769RSO/28XZ/LuAC4G9mAoCWLVg0PYuDYWs8m74B8JiT3HK/ny52LMHqZPBec2fmsB1BnNbD9aztFOmKn1NXKemH0ASZrVvORL4B7gBOYmyP5mo2PbjpKAOirxJZPJCA1SEKZWEqpWEqpWoJDJkLCVaZYGRbInRmtzyrkZl7uL9fdkGRzYsZlF08a4SHAYuJz4Yf+muOBXfJvZj0P6EO585AyakNNeFSw9UlJLTaOZGRfGn/MFSQWt8UsF8POZan441rZCjzf8khXKBysHIkmjkCvqkcnux2L+FzK5gvNTr0StCaHgUA71NVVeJwFBS1zBvCUrGRoTyfmDY4iPUBMbrkbVwbTaq+Za5s+fz9y5c3n99TWcOH2Gy4fHcKaqkflPvuz4ns4muzvLwFUxAMjkcpAsJI75jFv/OIF3n02g+HgIL/4xkSGjNrodQrijqLqRD7JOc+NF8USHBp31dRf0TfxOAXz9SyWVZu/Ecp4+Cw6LZdMbsez4NAaAYRfUo419iP3frUYZpMZiaiI6dmCrGH3nhyth+CiP03RhaiXffbqW4z/vI3/z29zx2GPExnpfbOP991ssjJdeetHxeRw2K6Go2sjhYluJ8Ibm6kTuLAPnIYNMJkeSrFx4eTqhEdHUVJYxbHQDf345jyWzDmBsmMWR7P8HKNi18QFHSvMz/81p9zpmmq1cNyaWIdGi/mAg4HcKoLTOhNLLOWr7eP2zf31KyamlFOZpkCskrr2zhEm3lvOfpbluU2k9pdg6L//xiw+x1lXy6M3jMBpbpsxef/11Xn/9dTQajU+ysWQyGQMiNQyI1HDliBhOGRo4eKaG3yx50THN6WwZuJPdeb0m1Mqf/xXGB8tf5dSRucBC4FJgJslXJ7d7He2Wwmc/FXH1ef24MCGiy99P4N/4XXvwlV8eRKluP269tWd9LvAiEAIcYeDwvzP3b3M77cEPCVIydmA4YwdGEKpWUlhYyMMPP+yothMcHMy1117L888/360FN6oaTOScquJQYQ2mDhx0roFEdXV1fPnGc+zedAz4EEgESoFZwHcOv0ZH+RHJg7VcMTy6WzIK/T3W3t/lA9/kAvj9LIA7Fr/1NRddMQuZ/ENgDRBCVNw2kq9+jMK8dZ3y4IcEKbnqvH7cddlgLk2MdmTQDRgwgIiIiFbVdsLCwrq92k5ksIqrz+vHbyYMQTc0CmU7barczVzUGipITRvEiOQFwFdAf2ALg85bxcNv2qoNd1SNOOekgY0Hi2kyixmCcxW/GwI4Y3+z3XTvYj59bZnjDVdVNphf9r2IZI0HaoA/UFm8jsrmIjhn06RDo1KQOkTLhQkRHp15JSUlzJs3j7lz57J27VoKCgp8+0U7kG9CUjQXJkSw53gl/yuscQwNvGlO8tFLT3I0ZwYy+TIk6yJOHbmXz/5VRFX5rdz5yNIOS47ll9XxcfYZrr8wnrAA6l8YKPj1L+pcYrvkZB6b31lF3JDn+XxNHBazjJCI49z2p8Ps21bKz7tVyORytzHw7pDLZFw0KJJLh0URpGzfEHJ24r3wwgu9UmknTK1kysj+JA+K5Puj5ZyoqPc4XTjl9gWO/WoNFUyYMYvLpg/n8zUvcnT/fH76Ph54jk9XreLoflt48lW33Omx5FhprZEPsk4z48J44iK87zwr6D4KyuvxRYK3XyoA9xV4tezeNBuwmd4Tb6jg+rsbUAUN4pcsLRazLSPQ05vMmSHRIVw5IqZPTndFhwZx08UD+KW4lu+PKNy+wZ3LhTlHOeYf/BWS9Q3gM2ACP+1IAPazb+t/mfPA0najIO3hw+mj+3NebBiC3sFildiZV8FPp6qYdUHXfwe/9AHYx6bKIPvb5lIgG7gFMDB01N8pzL+WJTOTHRFzdmwVe2VuY+DD1Eqmj43nposH9MmH35mRcWHccelgZI3VTPAyj8F2XQcBlwC7gCHAD0jS9SyaNoan7pjSbhSkPXx4j4dELEH3UtVgInPfaXJOGqgqL2XWrFkUFRV16Zh+aQG0zHebQPYXkJYCKmAPMJuC3OPNW8pImTTDbfCO8xtMhowLBoZzxfCYDs39voRGpeCbTZ9wsqKeLbmlDEwaBdgqFrnDfl2hGEjD1s7h18AG4CEk6TmgfR+KhMSP+RWU1TaRPqr/OXU9/Zlfimv59pdSR8WrzetXsWfPHpYuXcqLL77Ywd6e8atfr7CwkNV/nUd1RSmGUgva/nqQlgMqVOpXgSuA4057SGR/uxGTsdGj6R8ZrOLm5AFMGXnu3qyDo0P41fhBXpnmNp/AHC684krgTuCv2G6D5cjkawGlx1bo0JKhmP1LPh9ln6G6oW0xFoHvMJosfPlzMV//r5gmi9WpW9R7SJLE6tWru9QFyq+eiGXLlnH8UA7vPbeZI9mrMJSmEBxm4bePnyA17TvAhEzWIrJKrSGyXxyp6Te5NYEvTIjkV+MHMSjKv+rhdQcalYLrxsQxdXQsQe2EJtsbnkpWKxNmzOHPr+g4f9xKoBHJ+jtksq8xGYM9+lCcpxzLao28n3WaU5UN3fjNApeTlQ28u/dUq7Rt16nb4OBg5syZQ25ubqfO4YvWYDMBAzBOr9c/25ljaLVaGhsbAQXwN45k299K27n/1Xi0/c3s3Wx7c9VWVXBgx9fIZDaP/wWXTnLEudsj4kKClKSN6h+Q7bRHxYejGhPN3mIrpbWti344Bww5OwfnPT2KVQ/+hZOHn8JknIwm9ABlhYtb7dvelKPxvzlcPjyalMGiF6EvMFms7M6vZP/JKiSX5hPO4eBqtbrLXaC6ZAHYm4Lq9fotgKGzTUIPHTrETTfNA74FHsXWceMJYApP3zmSxTemtHlz/enlD906vZL6hXL7+EEB+fDbCdcomTluIGMHtg7lbS/V+ffP/oFFq0uIH9ZIY10Cp4++Qf7BFsupvaAhSZLYcbScTSJoqMsUVjXy3t5T5Jw0tHn47diGcbP59NNPmTdvXpe6QHXVApgNbG7+nAekA/vO9iADBgwgJkaJLZz3FPB/wPZmp951rebz3dXyB1DIZVwxIoaLEkTlWwClQs7kkf0ZqA0m9byEVv0MPTn5omLNLHjuOO88k0CuPpx/LR7KDfN/Yf/2O73qU3istJaKuiamj43r87MsPY3ZYmWXh7e+K3c99gIymYwLLghj0qRJXTpvVxWAFqhw+jumvY3bC6A5cyafW25Zx8GCPI7s2+4w8RUqDQp1iEfPNkCkRsmkYZHEqIw9EqTj722jnOWLksGXW77jwcee4Kdd3zge4NGXTuaaO//o9rrOWpTLV28NYc8XcXzy6mjgeja99RJ1VeWkXpOBbuot6DdvoLK0uM3+dXV1vPFdFZcNDef8/p59L33pGnY3Z6qM7DxeQ7XRfRs7d8hkMgwG77f3hN+0Bvvkk08ASJ00jQkz5rRp1+WJUfHhXH1evx738Pt7ooizfLGxsSSPSCB7e5MjNTosIpL4QUM97r9vyxBgPvA88DhZm9cDv0Opkpj9p78xfExKq+2d/Quh0f3JLjFTi8Tk8/uhVrkvMtKXrmF30GiysONoOYeKmkCpJlTpOcrSNeFLJpOh1YZ1WcauPjUGwF5QTgt0uZrmHQ8vJ2PhY47Cmc4mvzMKuYzJI/vbvN7n6PSeLykvK2PevHmseu9zrwqf2sb8+SiUM7HlW9yOXLGD3z35MdC2YKk7/8KRklre1Z9ut1VboJJbVMO6H09yqKim1XLX62rH9frWGBTs39/1YVZXLYD3AV3z5yRgSxeP5xXhGhXTxsSJuPSzwDmf4TJdClsOlWK2enbY2b3NVsuHKJSTsZg/wmrR8dbfyljwzyp2b7TdkE/9egqS03Fc/Qs1jbYy7qlDtFwyLAqF3H+blfYEpTVGvjtS5lEputZmcD/7EgSsZMCASPbvtxDRhbINXVIAer1+n85GOmDQ6/Vn7QA8W4bGhHDN6Fg0HsxKQcecFxtGhEbFxgNF1DV5HkfWGipAJsNizsIWPryBpsaJrPyDBmgAJHAJCXaXiCVJEvqCSvLK6kgf1Z+4CE23fC9/xmiysDu/kgNnqt33s/AwzapQqpyiXeORyV5HktIBGDGikbo6WZcUQJdtZ71ev1qv12/R6/Wru3qs9pAhY/ywKG64MF48/D4gLkLNbboEYsM9W1F3PfYCj/5nW/P0XzUwBZn8XSAM+BR4BGWQhpiBQ5DJZB0mYlXUNbF2834uuWIShUWdn7rqS1isEjknq3j7x5P8dLrKYw7FwpXvEhoZ7ch/UQapCY2M5o8vvEdQcBgm4++Bg7aHX1bBTffmsm5dKQMGdE2+PjF4DlLImTY2jssSu6c6TaASplaSkTKQEe2EELcuXQ6S9f8Ii1wOWIGnMDf9G4s5yOvGKpvXv8qBrB/50+P/4Fip96Xe+yJHS2p5Z89Jvj9aRqOp/QK3P37xAXVVFa2mWeuqKtj2/k5++v5x4J9AKOqQz0EaTcnJJ/HFo+CXyUDOaINRN6zWAAAT/0lEQVRVzLhQVKrtLlQKOdddEMvuYBX6gkq327jWITy4azmjdBJ5B/9MU+NtNNamMmmWleg4kyM2o013ZlcT94tMxgzORBWk5kRRKVEh587vW1Bez57jle12sLLjriybLWYjCFjM/u1/bf58Cvg9xvrPAdi18T2GDHmvy7Up/doCGBoTwm2pCeLh72ZkMhkTkqJJHxWL3M1rxR6FaZ+ZWfLOt8z9+/X88YUC+iUYaawfzgv3JZKrb5mudfVauzNxUybP4OE3v2L9nlN8d7iM2rOYB/dHCsrr+TDrNJ/9VOjVww9tIyyVQWrCo24A2T5s0bBB9Bu4id8v30fKJKlVJObNN9/c6RwAO36rAFIGa7nhwniPc8gC3zN6QDg3Jw9ArfTumscNaeI3j+sJDttBfY2SN5YM4aEZmSyaNrZNV6MX7pvVysS1mJocvgKrJPHT6Sr+s/skO46W09B09v0gegurVWLngaNcfNlVrPv2J68ffPt0H7T0gFCoYjA3PUdN5ScgjQGOApM4L+VtksaOaBOJ6YvalH6nABQyGemjYrliRIwY7/cCCdpgZqUORBus8mr7Hz59iYbaqxl8vm2a0WpdAmxEJmvtBJSklqlCT0VbzFYr2ScNvL37BN8fLafKj1ONm8xW9p0w8PaPJ3nkib9zeP/esypG62wh1RoqGJn6DMGhJ7CVcbeSMPwT7lt5kAkzBjiuk30odt/K9UyYMZuysrbl284Wv/MBTD1fy+gB4b0tRkATFRLErNQENh0s5rTBfaqv69j15OE5wFpgPTANScrCltPxPf0GDqGqvKRN3UJPkYhNFis5Jw3sP1lFYr8QLh4U6Rcp3ZIkccrQyM4DR/nbA7+nIPenDouyutJ2um878AK2ku0w6LwGZv3pDAOTzgdAG9tSpq11HswSbjsXS4LFhHr35hF0LxqVgpsuHsDoePfK2N3YNSr2ZxSqCdjKjQ0GviF2yH+oKCltt26hJyQk8srq2JBzhrd3n2BXXgXltW37GHQnkiRRXN3IrrwK/r37JJ/knOG1F54j72AWF191Xbtl1d3RUu4uBNvbPheYhUpt4Yb5Rdy3Mp+BSS2JW75uVuuK31kAAv9BIZeRPjoWbYiK3XmVrbLUXDsbW0xNBGmCsZqPoVBNxWJaDCym5MQdwFD6JfyVX//1zx4rD3fEydNnWLbQNqswbHACSf1CSdBqiI/Q+DwUvL7JQkFlIwcqSikor3cES7m+vbO2fub47E0xWrBdN1PTRZibHgfGAxAVt5ffPxtJVGyLE9Sbku++QCgAQYfohkahDVax2SV8uO304FbH38/fNwvJugVYB1xJ2emNrFxwPwrlxzzz3+x2szvd4RoiW1HXhL7AVt49NlxNfKQGbbCKqBAV0aFBXnU5Nlms1DSaqW40Y6g3UVzTSFGVkepGE3V1dYSGtvZBuCvDHhIeyYjky9otq95yvRR88VYsB3cuBSA8qp6BSa+iUn9FVGzrnBdPJd/bK3XfGYQCEHjFiNgwwl3Ch11rMzjXZ3j0P9v4/PXlHNh5KeamF4EMYA1DRz1LZUkZQR4SPDuMH3B5E1oliaLqxjbed6VcTpBSjrr5n0Iuw2KVMFkkLFYrRrPV0Yj1bDi6/8dWwxl3FalcsZhh18YovvpPLI11ChRKiasyykmbU4o6eDowvc0+rhaWN9ZFZ/A7H4Cg9yksLCQ9Pb1Nyem4CDW3pSbQP6zjJCz7DWwxnUGhuh34FUpVNXkHo1lx73D2fNkfq9PzZ58W2/TGP1uNeTtqX+YJs9VKfZOZyvomiqobOW1ooKi6kfI6I4YGU6ce/i3rV1FTWUbskCTuW7mecVNuZP/3X7XJ3LMjSXBwZzgr7h3Op68NoLFOwUhdLQ+sOsb035agDm6/8Iez1985wlKGjPN91JtBWACCNixbtoydO3e6LTkdplFy67iBbM0tbVWs0h2uQ4TK4rko1Ws5+EMEm9YMY/+3Ddz8hyKGjW5wZBXmH7Tt6/ym1029+azehK5WRFdp06im4BgrF9zaXKBWcgxLnDnxi4bP18SRf9Bm6iiDTpCxoBrdVLnXIbzuql8NigpmYlIMcRFqnxS/EQpA4KClOKuN1atXO8pOO4ebqhRyrhsTR7+woDbOQWfc3cCSdIoDP4Tzyap+nD4azCv3JwJvAv2A1jd0yuQZXH/3Ija88pTbdu6ecPUXdBXX8bgde2yDs7JauPJ/bHozksNZtlmOkAgzsYP+w/H/3cvJw7cw/prOyRMdEsTlI2J8XutSDAEEDg4dOsTs2bMJDrbNuXdUclo3NIoZF8a1W4bcFZkMLrqihhHJdwNPI5OZgN8Cx4C/A5HI5HJkMpnjTe8aiuypSExLzfzWUYiLb0xxu72n4huuuI7HAfoNHNJqWDJy/EKGX3yK5xcOb37464F/UF8dw/H//Q5ocsizaPpYxzk7kiFIIeeKETHM6aZCt0IBCBy4a4XeUcnpxH6hzEpNINrLZB77Q5q9bT3wKJJ0AbZehWHYKkLnEzt4LRdfdVu742t3D4672ITQyGjue/5dt8fwdo69uqKUn763NVG1ReHNwWKxYDIakStnYDJ+zC97X+IXfX9sdRL+ia0+zsNAteM4KnVz6nTzuduTQYaMUfHh3HHpYFIGa7utkIpQAIJW2Fuhb9++3euS09GhtsjB9tKK7bj2fVSpTxHZbz4jdQ8xeGQJEEVxwW848MMa6qv/ysY1H7o9jrsHx53nvK6qgt2bPmgtgwdL4eEbknl10V2cPpbLm0vuaVXurKG2iiC1hoFJo7j+7iUEh91PcFgBVvPnwDRk8kYuve4UYyfMQ6V+BChueeCbB/0mYyPlZ060Oqc7ayUmVE1GykCmjo4lVN29o3ThAxC0wrUVurcEKeVMGxNHToSGH46VY/VQ+KJldqDJZSrtLqCMh66/A6tlCRbzFcAi9n1jYd83G5DLV7Pss+d55Jb2pwXtVYzaK4N+8VXXkbX1M+RyBVarxTHHLpcr2Lftv4529K3LncnYtfEMuzbuxha2e7/t+8SYuPyGCi6dZsBiLuL5hV85FJDJ2EhVaRHj0m5k3JQb2PDy3ygvOoVktdqy/rQx1BjKHS3tL7o8nSf+vpS0lATkPVQ6TVgAAp+SPDiSW5IHEt7Om6vWUEHqNRlup9IeefsRUiY9g0J1JfAOtsIjM7Fav+bpO89jXNoxzh/3Z5RB7qcFW1cxar0NwKJpYxwRfNbmeUiTsZHsbzaStfUzW+jviWNIkoRklQNXAc9iy8zbDvwOCGfo6HrmLDrN4jePMGV2OaERljbThHFDRmA2NRGk1nB+ygRGJF8GktQqctJZEV40LJ6pqef32MMPPrIAdDrduJ6oByjoGwzUapgzfhDbfil1W/XnrsdeaI60CyVIo6GhpsrhsW+xEDYCO4BFwL3AHVRXJLHny0HYxtj3I5Nvx2TcCoQSHtUy3ecpiMbVmy+TKzg/dSJqTQgFh3KorarFYjofWzv6q4DrgCgnyU+SMFzP/z18If0HtVghnqYJ7dgtEJlM1iZy8qobf8X9C+5l00fvdLnVd2fwRW/AdGx9pod3XRzBuYJGpWD62HgOnK5mx9HyNhWIn/rVFR5N+VHjryI17SZqqyo4sm8nVuvjKIOWkXTh74mKvZ+srRrMTYOQrLcDt5P9DRzJNhI31EzsICOxg5s4k3ceF1y2hHFTruDnXRuoLCnCahmAJCVhMp5AoeyPxTyQ+qobiIiaisVSgsU0BGiddShXHCFuaC5X3BjJycOvUWsopf+g9sN23Zn3rm3rZ/1xCc88909Sh2hRKuRce9Ul3fI7dESXFYBer9+i0+nyfCGM4NzjwoQIErQaNh8qoaSmZVz+/179hG3vvOI21t3+kHz00pNIkrXZZDYSE3+GjIXVZCysprigimM/hbJ9wwkqi0dQa4im1qDm2H57jPFTAPy8C2wd6+DpOwFsgU2W5rybk4dt/+xvelXQCRSqfYREHGbIyDOYm352TDtecu2jbr+jp8QoZ/PeOXgpqV8oV4yIIdLLmgvdSY86Ab2JXPL3llHg/zL6o3xXD1ZxoNBE9ularJKEKjgMefPDoVQFuW0DZygr8diK7MX/52x2y4AhwGjk8rGkpC2iqjSIxnolxnoFjfUKTE1yVEFWlCoryiAr6mArkf2MRMU1cTL3MwoOfUTy5IHcvGABkEBjYwwajc2H4E3ikqusuXu+ayO70mLk0iHhDNLKMdZUUlLT4WHbP6cPfme/aQ3Wme16E3+X0R/li4+DcSOa2JJbQn5RBcba6jYRfs5t4OY++Yrjs70VmT3Md+HKd/ku8w0nC6KYsRPHcv3d1xMR7V26sevYPecbyPnm3yhVQTz67o52W9K54k5WOxdcrGP8MC0XJUT6fD6/q79zhwpAp9PNd7M4r7kluEBwVsSEBTFrXALfHTQT/eTLGM02T7ynTDpX7PP/P37xQZez5TpKue1qToFcJuOiQZGMH6r1214WHSqA7m74IQg85HIZY+JDuHRUDDvzysktrO2wJba7tGAAmUzOfSvXd6rQSHspt3V1dZ3OKZAh47y4UC5LjPaLcX57+GIWYKbtP91MvV6f6QOZBAFCSJCC9FGxjB0QwY5j5e02EW3vbR0R3d9rC8IV14zFmsqyTlfjkSFjRGwolwyL6jOl7H0xC5AJiAdf0GniIzXMHJdAQXk9u/MrWs0W2PGmQEZnTHZ3GYvVFaVsWLWMX/Z+51U1HplMxnmxoYwf2ncefDsiElDgNwyNCWG2bhDTx8YTF962gainAhl2PCXWeJv1Z8fbajxqpYLUIVp+c9kQrr0grs89/CByAQR+yPD+oQzvH0phVSP7T1VxtLQOSZLcvq2h4wKadsWw6Y1/UlF8xisLoc7N0MBOXLiG0QPCGRkX5vOCpD2NUAACv2VApIYBkRpqG80cKqrhcEktFXVty4J78g8c2LGZRdPGOLaz5wA8dccUnt10oN1zz3nwWcc0YMbCxwhTKxkZF8ao+PA++ab3hFAAAr8nTKNk/LAoxg+Lory2icMlteSX1VNR14SE1GHsf/a3G1sdT5KsLJo2pl2nngwZceEahsYEMywmlNjwoHOyU5VQAII+RUxYEBPCopmQFE1Dk4UzVY2cMjTwXq2BidfP4ZLrZjpMdrtiANt0ob2ElzunnlqpIDY8iLgIDXERahTGGoYmdK3vXl9AKABBnyU4SOHwF/zwtc28rzWaqcyYQk2jCaPZymdNNVw/+04MFeXs2LwRmVyOucnIsPh+3HblWCI0SsI1qjZ9BEpKzq5vQV9FKADBOUWYWkmYWok9q2/z5xsAmD17Nvfccw9z585l7dq1FBUVcZ6PSmv3ZYQCEAQEna10dK7Tt+cwBAJBlxAKQCAIYIQCEAgCGKEABIIARigAgSCAEQpAIAhghAIQCAIYoQAEggBGKACBIIARCkAgCGCEAhAIAhhfFAW1lw0frtfrH+rq8QQCQc/RJQuguS/glubS4UnNfwsEgj5CV4cASdgbr0Fe898CgaCP0KUhgEvTkHHA+562BdEbsKfwd/nA/2X0d/nAj3oD6nS6ccA+vV6/r73tRG/AnsPf5QP/l9Hf5QP/6Q2YLhyAAkHfo8u9AXU63Xy9Xv9s8+d00TRUIOg7+GIW4B86ne6YTqer9JFMAoGgh+iqE3ALEOUjWQQCQQ8jIgEFggBGKACBIIARCkAgCGCEAhAIAhihAASCAEYoAIEggBEKQCAIYIQCEAgCGKEABIIARigAgSCAEQpAIAhghAIQCAIYoQAEggBGKACBIIARCkAgCGCEAhAIAhihAASCAEYoAIEggPFFazB7Y5CpojKwQNC38EVR0FnNtQHHNfcHEAgEfQRfFAW1lwFP6qgxiEAg8C981RnoQeCejrZ7/PHHfXE6gUDgI2SSJPnkQDqd7kNgnl6v9/+magKBAOhiazD7mL/Z9M8D5gPP+lZEgUDQXXS1NVg6YB/3a4G9vhBKIBD0DF0aAuh0Oi1wW/OfqXq9vkM/gEAg8B985gMQ9A46nW4mYADG2Zu0etjuwfbWC/wfnU43ztNMm7f3gSs+mQXoLB0J3dkv1YPy2f0jw3sjCMrJB7NFp9MlebpBmuM1ptIL/hkvruE4IAlAr9dn9rB4dhm8vQ+TOuqW3V00/4b/Aoa7WefVfeCOXgsFdhYaMLgGEXW03g/kSwe2NN8QSU4RkT3JbGw3JticsL0hg0e8/A0XNz/4Sb0RSOblfZjXvD6vt4Ld7Of3sLrT90Fv5gJ0JHRv39wdnT/JaVle8989jRaocPo7xnWD5rfBFtflPUS717D5zboXQK/XP9tLgWTe3Gf/aP7fX4PdOrwPPNGbCqAjoTv9pXxEu+fX6/WrnczBcYC+pwQ7S6J78dwd/YbjgRidTjeuOZisN+jod96H7c1f6bLdOYHIBuwizSbhvl56MxhoecC1QLnzyl5++3tLuf3aNVsEfkXzTJcBWAa8rtPpesPS64h274P26E0F0JHQnf5SPsLb86f3Yhbk+7QMPZJozstovmnBNq6e2eysjO6F8WtH17CclnGtAZtF0NN0JON8YFmzc3Ae4DdKyul3dnsfeENvKoCObt5Ofykf0ZF86HS6+XavcW84AZ3enOmAwckK2dq8PtPJs651c4jupqNrmOm0vrcCyTr8ne00X8teCXVvto50LlaS/Xf2dB90SK/GATS/mfJwml7R6XRZer0+1dN6f5Gv+WJ/iG1cGE1LWrTACS9/4wpgfG9ZUl7I+GDz+ujemgbsLkQgkEAQwAgnoEAQwAgFIBAEMEIBCAQBjFAAAkEAIxSAQBDACAUgEAQwQgEIBAHM/wcVyBgkKNumfgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    # Initialize plot\n",
    "    f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "\n",
    "    # Get upper and lower confidence bounds\n",
    "    lower, upper = observed_pred.confidence_region()\n",
    "    # Plot training data as black stars\n",
    "    ax.plot(train_x.numpy(), train_y.numpy(), 'k*')\n",
    "    # Plot predictive means as blue line\n",
    "    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')\n",
    "    # Shade between the lower and upper confidence bounds\n",
    "    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)\n",
    "    ax.set_ylim([-3, 3])\n",
    "    ax.legend(['Observed Data', 'Mean', 'Confidence'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
